{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# lidisa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`lidisa` is a **li**ghtweight implementation of **di**rect **s**ampling in Python. This is a simplified version of the algorithm described in *The Direct Sampling method to perform multiple‐point geostatistical simulations* published by Mariethoz, Renard and Straubhaar in *Water Resources Research* (2010). It requires only `numpy` and `numba` as dependencies and makes use of the `numba` just-in-time compiler to significantly speed up sampling. `lidisa` is able to conduct both conditional and unconditional sampling of new images based on training data. This notebook shows a short workflow for using the direct sampler to create a new image based on a training image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from lidisa import dsampler\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "training_image = np.load('./data/sample_image.npy') # 40 x 40 pixel 2D Numpy array\n",
    "simulator = dsampler(training_image,sampling_mode='unconditional',radius=4,threshold=0.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`simulator` is a generator and we may use it to iteratively create new randomly sampled images. The state of the simulation is retained in `simulator`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "out = [x for x in simulator]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, many of the patterns from the training image are reproduced in the random realization!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAD6CAYAAAB57pTcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGmpJREFUeJzt3X/wZXVdx/HXmxUQhFhL1BACDUUatVAXbKxgzGghIVKSJkx3NFNKZf0xNtSIiDTajOn6Y5TCplVZFdTMxhFjTPMH/oDUTG3FUZcV/BUqLIv4A/X0xznfPHvv++6+z/187nnf7+7zMbMz3+/9cX7dc877e+557ftjTdMIAACMa7/sBQAAYF9EAQYAIAEFGACABBRgAAASUIABAEhAAQYAIAEFOIGZnWJmNy1o2jeY2aMXMe29hZltNrNLspcD+w4zW2Nmt5vZL9V87RzL8QIzu7T2dDGfu2QvwLIwsxsk3UvSTyTdLum9kp7RNM3tmcuFaWbWSLp/0zRfyl4W7J3MrH/cHyzph2rPDZL0tKZptgyZXtM0P5F0SO3XDtU0zYsXMV3MhyvgXZ3RNM0hkn5N0gmSLkhenn2OmfFHIdI1TXPIyj9JX1V3buj+TRVf9lvMgwLsaJrmm5L+TW0hliSZ2e+Z2afN7DYzu9HMLuo9d4yZNWb2JDP7qpl928z+uvf8Qd3XnreY2f9IWtefn5kdb2b/YWa3mtnnzezM3nObzey1ZnZV97XUNWZ2bzPb1E3vC2Z2wuQ6dK+5w8x+offYw8zsZjPb33n9QWb2hm6aW83s+f2vyc3sCDN7R/f+bWb2rN5zF5nZlWb2RjPb2a3Dwwe89+1mdrmZ3SZpg5mdaGYf67bHN8zsNWZ2QPf6D3Vv/Uy3Pc7pHn+Mmf1X956PmtlDevM4wcw+1S3bFZLuOvWhAwOY2SVmdoWZvcXMdkp6gpn9upl9vLffvmrlWDOzu3TniGO63y/vnr+q2y8/Zmb3Hfra7vnTzOyLZrbDzF7dnSM27Ga5N3c/H9vNZ4OZ3WRm3zWzp5rZSWb22W49Xtl77/3N7ANm9p3uHPcmMzus9/zDu2Nwp5m91czeZrueJ880s8900/2ImT2ozqexijVNw7+2HecNkh7d/XykpM9KemXv+VMkPVjtHy0PkfQtSWd1zx0jqZF0maSDJP2q2q+sju+ef6mkD0v6eUlHSfqcpJu65/aX9CVJfyXpAEmPkrRT0nHd85slfVvSw9QWjvdL2ibpiZLWSLpE0gdmrMd7JJ3Xe+4Vkl49Y/1fKumDku7erf9/95ZxP0mflHRht4z3k/QVSb/bPX+RpB9IOr1bppdI+viA994p6azutQd16/oItbdIjpG0VdLG3rI2ko7t/f5QSf8r6aRu/k/qtsOB3Ty3S3p2t63P7uZ3SfY+x7/V8a9/TPUeu0TSjySd0dtv13X74F26/fyLam9jqXuskXRM9/vl3XH98G6/vELS5XO89p7d+eL3u+ee0+3fG2asyyWSNnc/H9vN5zXdsXK6pO9Leqekw7vzwHckPbJ7/QMk/XZ3TN1T0jWSXtY9d6CkmyQ9o1uOP+yW46Lu+XVqz5nrumP0yZK+LOmA7M83dd/KXoBl+dcdZLd3O3Mj6d8lrd3N6zdJekX38zHde47sPX+tpD/qfv6KpPW95/5MPytuvynpm5L26z3/lt6Ou1nSZb3nnilpa+/3B0u6dWI9VgrwOZKu6X5e083nxBnr8/9Fsfv9T3vLeJKkr068/gJJ/9T9fJGk9/We+xVJ3x/w3g/t4bPZKOmdvd8nC/DrJL144j3XSzpZ0m9J+rok6z33UVGA+Rf8p9kF+P17eN/zJL2t+9krqpf2XnumpM/N8donS/pw7zmT9A0NK8D36j2/Q9Ljer+/S90fEc60zpZ0Xffzo5zj/OO989hlkl448fyX1RX3ffUf9y12dVbTNO8zs5MlvVnSPSTdKklmdpLaq8QHqf0L8EBJb5t4/zd7P9+hnwUpjpB0Y++57b2fj5B0Y9M0P514/j6937/V+/n7zu+zAhvvknSpmd1P7V+vO5qmuXbGayeXsf/z0ZKOMLNbe4+tUXtVv2Jy3e9q7X2xyHv785KZPUDSy9X+xX+w2hPSJ2cs98ryPcnMntl77IBunRpJX2u6I77T3/7AvCb32wdK+ju13+Cs7Lef2M37Z50vhrx2l+O2aZrGBv4Pi6ZpQucXM7u3pFdJeqSkQ9Ve+d/cW47J+U6eQ841s2f3HjtAu57n9jncA3Y0TfNBtVeeL+s9/GZJ/yrpqKZpDpN0qdq/NiO+ofar5xX9/17wdUlHmdl+E89/beBiT2ma5geSrpR0rqQ/kfSmPSzjkb3f+8t7o6RtTdOs7f07tGma0wOLEXnv5JBcr5P0BbVJ559T+/X87rb1jZL+ZmIeBzdN85Zuve5jZv33V//vHdgnTe63f6/29tKx3X57oeLniHntctx2+/miitrfqr219uBu/TboZ+s3ef6Qps8hL3KO0SsXtKyrAgV4tk2SfsfMVoJYh0r6btM0PzCzEyX98YBpXSnpAjO7u5kdqfZr5BWfkPQ9Sc83s/3N7BS195XeWrwGrTeqPVDOVPtVVmQZ76P2Xs6KayXdZmZ/aW1Ya42ZPcjM1vmT2sU87z1U0m2Sbu+uKs6beP5bau+xrbhM0tO78IiZ2d2sDc0dKuljkn4s6VlduOWxkk4MLDcw1KFqv8L9npkdL+lpI8zz3ZIeamZndN84na/2/u0iHKr2XLXDzI5S+xX7io9IWmNm53XH2ePUfhOw4h8k/YWZreuO0UO6Zb7bgpZ1VaAAz9A0zc1qi9cLuof+XNLFXeLxQrUFK+pFar/23CbpavWuRJum+ZHa4nia2qDFayU9sWmaL5SuQzf9ayT9VNKnmqa5YTcvvVjtV0jbJL1P0tvV/rWrpv1/iWeoTYVv65bz9ZIOc6e06/znee/z1P6Bs1Ntcb1i4vmLJL2hS1M+vmma/5T0VLVhklvUhto2dPP/kaTHdr/fova++D/vabmBOTxXbQBwp9qr4cn9trru6+Nz1N6y+Y6kX5b0aXXHbmUvVPvH6w613wa+o7ccP5T0B5KervY4e7zaEOjKOeQTav+Qfl33/BclPWEBy7iq2K63xrA3MrP3S3pz0zSvH/Ce89SGyE5e3JIBqMnM1qi9rXV20zQf3tPrF7wsn5S0qWma3d362qdxBbyX677qfaj28Ne4mf2imT3SzPYzs+PU/jX/zjGWEcD8zGy9mR1mZgeq/cbux2pv/Yy9HKeY2b26r6CfIumBar/xwwykoPdiZvYGtf+/9vymaXbu4eUHqP3a7L5qk99vVft1OIDl9huStqg9hj+v9n9zLOIr6D05Xu0f+ndT+1+MHjeRsMYEvoIGACABX0EDAJCAAgwAQIJR7wGf+N4L+L5b0p3/Mv3f9PY/62bnlfNP79bjpzf14cd9e+55jKX2timdxxjLc+36lyy6WUN1Jzz95SnHcubnNIabr7/H1GPR49Z77xAl8yk5t0Q/u9XwGXvL+OlLnzPz+OYKGACABBRgAAASUIABAEhAAQYAIMGoISzvBrVnrABRNLRQO5xQO3DlTq8wkBEVXefwtg6GLzzR7RoNeGD5zPqcli2MM6+SY6d0PllhzpLjtlRJsMs99w1cRq6AAQBIQAEGACABBRgAgAQUYAAAEqSPhuTd8I7GYUoDQLUDD2u3Og1Pjpt7cq5ogKgkzDRrPp7oOpcEN6LLMlannNqhMMR5oSApfs5YjRYRzHLfXxC4qh0U884rmcdT9Bw7FFfAAAAkoAADAJCAAgwAQAIKMAAACUYNYbldVQqmF+3m4oULpPrBgZIOKp7w+gV5298NUQ1QEowo2f61O4wNeW/t4diAZRANe3nnjGggqSTAuIjwY3ZXPK6AAQBIQAEGACABBRgAgAQUYAAAEqR3wqodaPFuqq+dcVN9VledSWOEbqLTK+kK5QYlKnfqkuLbq/Y2HCscVTJNN+CxvmBh9jFDtv1YndEyeNvh3KOvK5rmlu3r5n5vSXCp9nCEYw1ZGV7n3RzfXAEDAJCAAgwAQAIKMAAACSjAAAAkSA9h1TYkeOHdQi8d4jDy3o2nXjX1WEkAwhPt8DJr25SsX1ZXqOrhKJUFN7xpRoN/mNG5bcDntEyBq5Kub975YizRTnnR/Trr3DfW+4d+zlwBAwCQgAIMAEACCjAAAAkowAAAJBg1hDVGOGfI8FJZgZjaoYNFGOOzWqZORaXzJXBVnxt+lB9yiQYqS+Z92MUHTz2248I75p5Hic2vPH3qsQ3nv6domm4nrfNj8/bCVZuuPi30XneY2lUwvGeNzntcAQMAkIACDABAAgowAAAJKMAAACSgAAMAkGDVtKJcxFivtVtMlqQuo0qSw97rZi1z7RRiNCU8/6iivrHGCHY/gxH2h73ZGMeTFN8foonn2su9iMRz1LufcvLUY/fQ96Ye8xLPUdHtP9axHD3Hem0679w6bLxvroABAEhAAQYAIAEFGACABBRgAAASrJoQVmZrspJQxVjBgSzRwIL3WEngaqyATomsoN6yi7aL9UIuQ1qGjhHuGePzzGxp+ph//ODUY14obO3W2PTcdQlu/7HOmyWh1qG4AgYAIAEFGACABBRgAAASUIABAEiQHsKKBiDGCjOVTDMayIi+zhtjc8tZ848l7AVdDh8QJKg9fm+4O1ZBIGPZQk97UwBvXiVjJ9fulCYt12eybME9L3A1pKPeJC9YF/3sS7oUznp/dhcuroABAEhAAQYAIAEFGACABBRgAAASpIewojeyFxGUqH1jvXaAa8v2+QNXpbygTG1jDOFXGmop+Uxrh9ZWo7G6XnnGCG56QUlPdLi+6L7pTS+6LEOmqWDnqigvcBUd1m9tdCbH+Q8vU9huBVfAAAAkoAADAJCAAgwAQAIKMAAACdJDWJndSWrflC/p6lV7vqWhlpLORMu0HaIWEdAIdwxaX33WKUqGcVtE6G+Mrnabrp8OLkU7t5V0xBvimXffPvXYq285euqx6Dmj9nFbe7jFtTP2pej+OeYQslwBAwCQgAIMAEACCjAAAAkowAAAJBg1hFX75n1mZ5NlC/zUFl1GNzwTDXCtgu1Q2764zvMq7SQWfX/tY7n2eS7aRWvWvuWFsLxpekOTLttwnpPcdZ7RCSsa9ItuB4YjBABglaIAAwCQgAIMAEACCjAAAAlGDWF5nVZcM26iz2tWcKMk5DFGmKbkxn+021DpkG/e+2v3NCrZDqUhmZIgz94cuCrZb9zObZXnIfkdlrx9s6Tz3mrgdb3ylAQqPd75vnbXqyGi+5M7XKbzuptVfsxzBQwAQAIKMAAACSjAAAAkoAADAJBg1BBW7eCGx7sJvoihzhYRaJoUXpdgsGFIF6BlChCVbIfofjNrH3FDI6s0jFNT9QBjMHg5VhAqOp/o62oHkoYcn14nrC3eMbV1+jOtHVCLbofo6xYRfox+Lt4yesGs3Q03yhUwAAAJKMAAACSgAAMAkIACDABAgvThCEtumJd0hSpVMqxZyTpHAwvRUMusZSn5rJZqOzjc/WHW9vK6NiV1dNtbeAEut2uVs38M2WeyhgyNrt9YvE5Y0W5knuh2LQkzucdo8Jw0K1BZu4OdF1obiitgAAASUIABAEhAAQYAIAEFGACABKOGsDzRm+AbT71q6rEt29dVX57aHX5qd+4p6SI0lnCIoWBbl3QWGhKdqL293f1hN51yllVJsDEasClV+3iM7tdjBK4W0QEqOoTp4QVhJu91tYeNHLJv1u74Fw5odrgCBgAgAQUYAIAEFGAAABJQgAEASJAewoqKBq6GhBPcQEYwEBBVu6POIrp6RS0i+FHTIpalpGvTvqYkiLOIeUeNtYwR0WUZsr9tuvq0uZenpOPfGOeLIefIkg597roEQ2u7wxUwAAAJKMAAACSgAAMAkIACDABAglFDWCXdixYRiqgdCKg9hJ93k3+M+Q5REhSLdt4pCZ5FpzdkCDOva1NJRzfpue68l1nWUH9jiQ5zN0YosnR71d7eJUGxMfab0nNaScgsOqTjCq6AAQBIQAEGACABBRgAgAQUYAAAEowawvICV14w686t02EH7+b2WF2hand0yeqysxqMEeCKdkCTVH2YPK+j2/nHV53FKMYIUpUedyUhoKzAVck+PKTj3xhDhnq8833t7Tqo82FlQ4ei5AoYAIAEFGAAABJQgAEASEABBgAgQfpwhENvWvetHSkoEQ0slLzO7cR0XGhyCwnEhLvOyOkGU/lziYZIMj/7fS1YN0ZIKXycSG6IbpkCV+HQkxNALZne0NdGLNO+PqSr3cyg5YQxh6fkChgAgAQUYAAAElCAAQBIQAEGACABBRgAgASjpqBrp8vCbQYLect97tHXTT226erTQu+dmdKryB2r0mkDN2vM4ZLt6Cbbk5KTtduIrpZ5L1rJsVeyDWb9rwlv346KJp5rr1/WGOCl74+ex0vWL2pIgt0765Z8pjVqDVfAAAAkoAADAJCAAgwAQAIKMAAACdJbUXpKbniPFUTwAldR0dBHNKwVnZ4XYNnohMkkadP1sUBZdBtGQ2GekrF/F6EkfOF+pusLFmaJjNHCb+bxGWzb6hkjkFRirGDpGGqf04a0oixpexw1dLtyBQwAQAIKMAAACSjAAAAkoAADAJBgKUNYnswxb2u/NyoaMHBDB8HAwZAwWe0wSDQUMda4z6iv5DhZRNAuug/Xfl1Jd6xwl6+CINos0XNQuMOVs4ze9BbRRS4tMLebkCVXwAAAJKAAAwCQgAIMAEACCjAAAAmWMoQVDV9Eb9TP6oxSMs2SkEBJuMQLLkVDGkNCDN4yRkNTJcEGb11qB67GCni481bOsIyr0SKCdiWf8xjDS7rTKwxXlQTFxhg61ZM5bGd0e5UMgbmCK2AAABJQgAEASEABBgAgAQUYAIAEo4awwiGGYOggsztWZkhgkheOKum8I/mdbbyOVCVDfEU/+2ggzJveIrrfLHsHtTHUDiRFO0DN2t9qzztLyb4+RPT9WUOnjqVoO1boPMYVMAAACSjAAAAkoAADAJCAAgwAQIJRQ1i1O0pFDbnxX7I8aZ1yRppeSeDKE+7Q43x+G4++buqxIUMrlhjjc152tcNt0Q5oOX2ZunlX7nQX7ojnbNch++DGU6+aemzL9nXua+dVMnRqSTArq1NXLVwBAwCQgAIMAEACCjAAAAkowAAAJEjvhJWp9vKUBHFqhxM8meGh2gG8zMBVCTc0sr7qLJZKeP+q0FVoUnT/qh2uig5zKi/UuIBzZO3AVVR06NSSc180/LWsuAIGACABBRgAgAQUYAAAElCAAQBIkN4JawyzgjReIMBVOSBSEq6KBktWQ7emrGH9vG0z1j6ybMOx1VQ7zDQkAFcSrvJEgzy1j7PSYURrqx0OrR3MKg1cZQ+jyBUwAAAJKMAAACSgAAMAkIACDABAglFDWNWH6QqaGWxYQPedSbXXORrSKA1z1A5slQRvxgiZZe4jq1HR0IPB140VGvT2TS8s5Aa9nOllBSDHOpajwucgTS93NJjlGVIrop/zonAFDABAAgowAAAJKMAAACSgAAMAkGDUEJZ7w9t53bJ1CyoJVdQeVqv2tikNHIS7H22NBShqh8xK37saOopliAbjvDBNSfiuVDRwFTVGR6rSeYwx3F/t81K0Y5ZnyPZypxkMXobns5vhRrkCBgAgAQUYAIAEFGAAABJQgAEASDBqCCuqJHTgBTcOu/hg97U7Lrxj7vmUhEaiIYZosKEk/DJESXetwwvWpWi+Cwjy7E1DQs7L24ejHa4WoahLUvB8M9YQgMsueh6pHXiLBrPCQ4iq7BzrGTJviStgAABSUIABAEhAAQYAIAEFGACABEsZwoqKBlqGhK2iN9trB65qG6ubWEkwpXaoZdnCWvuaknBONKQ3K2xVe7jS2mHA0YbP9JbHG+4v6bxUIhpwGhL0im6H6D4yFFfAAAAkoAADAJCAAgwAQAIKMAAACZYyhFUSdhhi46lXTT226frTph6LDlkV7cYTFQ1ueOuxZfu6uac3RMlnFe4aExweLDrfktdJBLYkP+iy8ejrph7bdLVzPAXnUft4GqL2Z1wyfKMXUCs997lBpcodCKNdx8LnEC9MNrDz1KTodvCW0f1MnWXcHa6AAQBIQAEGACABBRgAgAQUYAAAElCAAQBIMGoKunY7r9K0tJvQLEjVlbTiK2kD5yWePZnpXXfeTrrZTWoH51EyTu8ihOe9foSFGUHW8SSVHVO1j8eS4yycLnf+t8as+ZYcF7X/p0JJC9/S8aaXqZXoCq6AAQBIQAEGACABBRgAgAQUYAAAEqyaVpQlLQUH3SwfocVk1ribpe0Wa4+RWvu90ZacQ/avsdqirjbhY2qE40kqO6aWfRxcL9zmWcQ+GB1bd4ymodHPadayePudt35e8Mx7XVForcMVMAAACSjAAAAkoAADAJCAAgwAQIKlDGF5N8u9MTE9e/tYrSVhpkEBASc8s+zbNtoRbBGWfdvsi6Khm9qfXckxWtIpKhqYkuoHthbVKWoes0J+0e3jhb3WBrf3kM9A4goYAIAUFGAAABJQgAEASEABBgAgQXoIq2S4sugwYqUdoErmHeUGB4I39KPLMiQg4IUOoutXO4QyRpepIUO5lUxzbxmOMLpdso6nme93lrs00FRT9NiJrlupaFeokuO29j6yiM/Om48X9Rp6vuAKGACABBRgAAASUIABAEhAAQYAIMGoISz3BnVwOKgSQzqyLCIMEhG9ye9yulZFzdw2BcPIeV3LanfKGWN6s2QFxZZJZriqhLfcd24N7tfBz909fwU7y3nT8x778jmXTj32CJ09PZMZ74+KDteX1TmvJLwqlQ036m2HteE5t7gCBgAgAQUYAIAEFGAAABJQgAEASGBNk9PxBcDqdeJ7L5g6cURDepkhrKyuVyVhn9J5RD+XcOAq+N6sENasz9Nbl5LuflHbzn/uzI3IFTAAAAkowAAAJKAAAwCQgAIMAEACQlgAACTgChgAgAQUYAAAElCAAQBIQAEGACABBRgAgAQUYAAAElCAAQBIQAEGACABBRgAgAQUYAAAElCAAQBIQAEGACABBRgAgAQUYAAAElCAAQBIQAEGACABBRgAgAQUYAAAElCAAQBIQAEGACABBRgAgAQUYAAAElCAAQBI8H8km26FKe9htAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,axes = plt.subplots(1,2,figsize = (8,4))\n",
    "axes[0].imshow(out[-1],vmin=0,vmax=6), axes[1].imshow(training_image,vmin=0,vmax=6)\n",
    "axes[0].set_title('Randomly generated'),axes[1].set_title('Training image')\n",
    "[ax.axis('off') for ax in axes];"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:geopandas_p36]",
   "language": "python",
   "name": "conda-env-geopandas_p36-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
